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1 . 0  INTRODUCTION 


This  report  is  a  summary  of  the  work  carried  out  by  the 
Scientific  Studies  Corporation  (SSC)  team  in  Phase  I  of  the  "Two- 
Dimensional  Processing  for  Radar  Systems"  Small  Business 
Innovation  Research  (SBIR)  program  for  Rome  Laboratory  (RL)  . 
Specifically,  the  formulations,  analyses,  and  simulation  results 
for  space-time  processes  obtained  in  the  program  are  presented. 
The  report  includes  the  work  carried  out  by  the  University  of 
Central  Florida  (UCF)  under  contract  to  SSC. 

The  model-based  multichannel  detection  formulation  is 
presented  in  the  context  of  a  two-dimensional  (2-D)  representation 
for  space- time  processes  in  general,  and  airborne  surveillance 
phased  array  radar  systems  in  particular.  For  the  phased  array 
space- time  adaptive  processing  (STAP)  problem,  such  a  formulation 
requires  2-D  parametric  models  for  the  channel  output  process 
under  each  one  of  two  mutually-exclusive  hypotheses.  Results 
presented  herein  demonstrate  the  applicability  of  2-D  model 
identification  algorithms  and  methods  to  the  phased  array  STAP 
problem.  Specifically,  it  is  demonstrated  that  2-D  model 
identification  algorithms  provide  representative  models  for 
airborne  surveillance  phased  array  radar  simulated  data. 
Furthermore,  the  identified  models  can  be  used  to  design 
hypothesis  filters  for  use  in  the  innovations-based  detection 
methodology  resulting  from  the  binary  hypothesis  testing 
formulation  for  moving  target  detection  pioneered  by  Metford  and 
Haykin  (1985)  and  extended  by  Michels  (1991)  to  the  complex-valued 
multichannel  case.  This  is  a  novel  feature  of  the  work  reported 
herein. 

The  2-D,  least-square,  frequency  domain  (2D-LS-FD)  technique 
formulated  by  Mikhael  and  Yu  (1994)  was  adopted  in  Phase  I  as  the 
baseline  2-D  model  identification  method.  This  technique 
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approximates  a  given  2-D  complex-valued  field  with  a  2-D  rational 
function  model  of  the  auto-regressive  moving- average  (ARMA)  class. 
Of  course,  the  ARMA  class  includes  the  auto-regressive  (AR)  and 
the  moving-average  (MA)  classes.  In  this  technique  the  ARMA 
coefficients  are  obtained  as  the  solution  to  a  set  of  linear 
equations.  Excellent  results  have  been  obtained  in  the  image 
noise  canceling  problem  (Mikhael  and  Yu,  1994),  which  is  related 
to  clutter  cancellation  in  radar  space-time  processing.  This 
success  is  duplicated  herein  for  the  STAP  problem. 

The  2D-LS-FD  algorithm  co-developer,  Prof.  W.  B.  Mikhael, 
Chairman  of  the  Electrical  and  Computer  Engineering  Department  at 
UCF,  provided  support  to  SSC  during  Phase  I.  Dr.  Mikhael  is  a 
Fellow  of  the  IEEE,  and  has  made  extensive  contributions  in  signal 
and  image  processing.  In  Phase  I  Dr.  Mikhael  supervised  the  work 
of  Dr.  Q.  Z.  Zhang,  a  post-doctoral  appointee  at  UCF  with 
expertise  in  STAP.  Drs .  Mikhael  and  Zhang  are  both  scheduled  to 
continue  to  provide  support  in  the  proposed  Phase  II. 

This  report  focuses  on  the  phased  array  STAP  problem  for 
airborne  surveillance  radar  systems  since  such  is  the  main 
application  of  interest  at  RL .  However,  other  contexts  are 
possible,  in  radar  as  well  as  in  other  applications.  As  an 
example,  the  complex-valued  amplitude  of  the  signal  collected  by  a 
radar  system  as  a  function  of  range  and  azimuth  can  be  viewed  as  a 
2-D  image;  the  appropriate  context  for  such  a  case  is  then  space- 
space  processing.  Two-dimensional  modeling  and  identification 
methods  can  be  applied  in  a  wide  variety  of  military,  government 
non-military,  and  commercial  areas.  Specifically,  such  areas 
include:  radar  array  surveillance  systems,  active  sonar  array 
systems,  optical  sensor  systems,  non-destructive  inspection  (NDI) 
systems,  communication  systems,  speech  processing,  medical 
technology  imaging,  and  geophysical  tomography. 
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An  important  distinction  in  the  context  of  radar  system 
applications  is  that  the  random  processes  are  represented  as 
complex-valued  processes  in  most  cases.  Many  time  series 
techniques  and  models  have  been  formulated  for  complex-  as  well  as 
real-valued  processes.  However,  the  majority  of  the  open- 
literature  results  for  2-D  algorithms  have  been  established  for 
real-valued  processes.  This  is  true  in  particular  for  the  2D-LS- 
FD  algorithm  selected  as  the  baseline  in  Phase  I.  Thus,  the  2D- 
LS-FD  algorithm  was  extended  to  handle  complex-valued  processes . 

1 . 1  Problem  Statement 

Consider  a  linear  array  radar  consisting  of  J  equally- spaced, 
identical  antenna  elements  (or  identical  beamformed  subarrays)  in 
a  side-looking  configuration  on  an  airborne  surveillance  platform 
moving  at  a  constant  speed  in  level  flight.  The  array  is  aligned 
with  the  aircraft's  longitudinal  axis,  and  the  aircraft  velocity 
makes  a  crab  angle  y  with  the  aircraft's  longitudinal  axis.  The 

radar  array  is  radiating  a  coherent  pulse  train  at  a  constant 
radiation  frequency,  and  at  a  constant  pulse  repetition  frequency 
(PRF) .  Each  antenna  element  (or  group  of  elements)  is  referred  to 
as  a  channel,  and  the  ith  channel  output  (after  pulse  compression, 
demodulation,  and  sampling)  corresponding  to  a  single  range 
resolution  cell  is  a  complex-valued  discrete-time  sequence  denoted 
as  {Xj(n)}.  The  J  scalar  sequences  are  concatenated  to  form  a  vector 

sequence  {x(n)}.  Process  {x(n)}  is  assumed  to  be  stationary,  ergodic, 
zero-mean,  and  Gaussian-distributed. 

The  received  signal  in  an  airborne  surveillance  phased  array 
radar  can  be  referred  to  as  a  space-time  process,  wherein  the 
spatial  connotation  arises  from  the  spatial  diversity  of  the 
antenna  array  elements.  In  general,  the  received  signal  contains 
a  target  component,  as  well  as  receiver  (broadband)  noise,  jammer 
noise  (broadband  interference) ,  and  ground  clutter  (narrowband 
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interference)  components.  The  one-dimensional  (1-D)  multichannel 
(vector)  representation  of  the  radar  space-time  signal  has  been 
adopted  extensively  in  the  development  and  analysis  of  optimum 
joint-domain  adaptive  algorithms  and  sub-optimum  block-covariance 
algorithms  for  STAP,  which  encompasses  target  detection  and 
interference  rejection  in  airborne  surveillance  radar  arrays 
(Brennan  and  Reed  [1973];  Jaffer  et  al .  [1991];  Ward  [1994]). 
These  algorithms  are  often  referred  to  as  the  conventional  (or 
classical)  STAP  algorithms.  More  recently,  Michels  (1991)  and 
Roman  and  Davis  (1993a;  1993b)  have  adopted  the  1-D  vector 
representation  for  the  development  of  multichannel  AR  and  state- 
space  models,  respectively,  for  joint-domain  innovations -based 
detection  algorithms.  Alternatively,  a  space-time  signal  can  be 
represented  as  a  2-D  scalar  process.  The  2-D  representation  is 
essential  for  visualizing  and  understanding  the  spectral  energy 
and  correlation  characteristics  of  the  total  signal,  and  of  the 
ground  clutter  component  in  particular.  This  representation  also 
offers  algorithmic  and  modeling  advantages. 

With  respect  to  model-based  detection  for  the  airborne 
surveillance  phased  array  radar  STAP  problem,  1-D  models  offer 
dynamic  and  static  degrees  of  freedom  in  the  temporal  axis,  but 
only  static  degrees  of  freedom  in  the  spatial  axis  (see  Sections 
A. 2  and  A. 4).  In  contrast,  2-D  models  offer  dynamic  and  static 
degrees  of  freedom  in  both  axes  (space  and  time) .  Thus,  a  model- 
based  detection  methodology  using  2-D  scalar  models  is  inherently 
better-suited  to  the  cases  wherein  channel-to-channel  correlation 
exhibits  a  complicated  structure.  Such  cases  occur  physically 
when  the  platform  and  scenario  parameters  lead  to  a  non-integer- 
valued  clutter  ridge  slope  parameter,  to  non-zero  misalignment 
between  the  array  longitudinal  axis  and  the  platform  velocity 
vector,  and  to  non- zero  internal  clutter  motion.  These  conditions 
are  common  in  airborne  surveillance  radar  scenarios. 
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1 . 2  Report  Overview 


Model-based  multichannel  detection  in  the  context  of  the  2-D 
representation  of  the  channel  output  vector  process  is  discussed 
in  Section  2.0,  including  a  summary  of  2-D  time  series  models. 
The  2D-LS-FD  model  identification  method  and  its  formulation  for 
STAP  are  discussed  in  Section  3.0.  The  airborne  surveillance 
phased  array  radar  application  is  discussed  in  Section  4.0, 
including  simulation-based  results  generated  using  the  SSC- 
developed  radar  data  generation  software  described  in  Vol.  II  of 
(Roman  and  Davis,  1997).  A  summary  and  suggestions  for  further 
work  are  presented  in  Section  5.0.  Appendix  A  presents  a  2-D 
representation  for  1-D  vector  state-space  models,  which  allows 
comparison  of  results  obtained  using  state-space  model 
identification  algorithms  (Roman  and  Davis,  1997)  with  those 
obtained  using  the  2D-LS-FD  algorithm. 
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2 . 0  MODEL-BASED  MULTICHANNEL  DETECTION 


The  model-based  approach  to  multichannel  detection  involves 
the  generation  of  an  analytic  model  for  the  multichannel  system, 
utilization  of  such  model  to  process  the  multichannel  data,  and 
determination  of  a  detection  decision  by  applying  a  detection  rule 
to  the  processed  data.  A  model  class  is  selected  to  represent  the 
multichannel  system,  and  the  model  parameter  values  determine  the 
specific  system  model.  Model  parameters  can  be  identified  on¬ 
line,  by  applying  an  identification  algorithm  to  the  received 
channel  data.  Alternatively,  model  parameters  can  be  identified 
off-line  for  various  conditions  and  stored  in  the  processor  memory 
to  be  accessed  in  real-time  as  required. 

For  the  work  reported  herein  the  time  series  class  of  2-D, 
linear,  shift-invariant  parametric  models  is  selected.  This  class 
includes  MA,  AR,  and  ARMA  models.  Emphasis  is  placed  on  the  ARMA 
model  subclass  because  it  is  more  general  than  the  MA  and  AR 
subclasses.  One-dimensional,  vector  (multichannel)  models  have 
been  applied  to  the  multichannel  detection  problem  by  Michels 
(1991)  and  Roman  and  Davis  (1993a;  1993b),  and  the  success 
obtained  in  their  work  provided  motivation  for  this  research. 
Utilization  of  the  2-D  model  class  in  the  context  of  multichannel 
detection,  however,  is  the  most  important  and  novel  aspect  of  the 
work  reported  herein. 

Two  types  of  2-D  model  parameter  estimation  algorithms  can  be 
defined:  (a)  algorithms  which  operate  on  the  channel  output 
correlation  sequence,  and  (b)  algorithms  which  operate  on  the 
channel  output  data  directly,  without  calculating  correlation 
information.  The  AR  model  identification  algorithms  discussed  by 
Therrien  (1986)  are  examples  of  the  former,  and  the  2D-LS-FD 
algorithm  formulated  by  Mikhael  and  Yu  (1994)  is  an  example  of  the 
latter.  The  2D-LS-FD  is  the  algorithm  selected  for  this  program. 
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As  the  name  implies,  this  algorithm  is  based  on  the  least-square 
philosophy,  in  contrast  to  the  stochastic  philosophy  prevalent  in 
most  of  the  model  identification  literature. 

2 . 1  Multichannel  Detection 


Detection  problems  in  the  context  of  radar  systems  (as  well 
as  in  other  applications)  can  be  postulated  as  a  binary  hypothesis 
testing  problem.  Specifically,  the  class  of  problems  addressed  in 
this  report  involve  the  following  two  hypotheses : 

H0:  Target  (desired  signal  component)  is  absent 

H1 :  Target  (desired  signal  component)  is  present 

H0  is  referred  to  as  the  null  hypothesis,  and  H1  is  the  alternative 
hypothesis.  The  model-based  approach  to  detection  with  phased 
array  radars  is  couched  on  the  assumption  that  the  process  at  the 
array  output  can  be  described  by  a  linear  parametric  model  under 
each  of  the  two  hypotheses,  and  that  the  parametric  models  that 
correspond  to  the  two  hypotheses  are  sufficiently  different  to 
allow  selection  of  the  correct  hypothesis  by  the  evaluation  of 
measures  that  are  sensitive  to  those  differences.  In  innovations- 
based  detection,  as  formulated  by  Medford  and  Haykin  (1985), 
Michels  (1991) ,  and  others,  two  parallel  paths  are  established  for 
processing  the  received  radar  signal.  This  is  illustrated  in 
Figure  2-1  for  a  general  case  with  off-line  parameter  estimation. 
In  the  top  path  the  received  radar  sequence,  {x(n)},  is  processed 
with  a  filter  tuned  to  whiten  the  sequence  when  H0  is  true,  and  in 

the  bottom  path  the  received  radar  sequence  is  processed  with  a 
filter  tuned  to  whiten  the  sequence  when  H-)  is  true.  The  filter 
outputs,  {v(nlHj)  I  i  =  0,  1;  n  =  0,  1,  .  .  .  ,  N-1 } ,  are  labeled  as  residuals; 
these  residuals  are  true  innovations  processes  when  the  data 
condition  (hypothesis)  matches  the  filter  type.  This  architecture 
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is  referred  to  herein  as  the  multichannel  innovations-based 
detection  algorithm  (MIBDA) . 


{x(n)} 


Figure  2-1.  Multichannel  innovations-based  detection  architecture 
with  off-line  hypothesis  filter  design. 


An  alternative  to  the  MIBDA  detection  approach  is  the  model- 
based  or  parametric  adaptive  matched  filter  ( PAMF )  formulated  by 
Rangaswamy  and  Michels  (1997).  In  essence,  the  PAMF  consists  of 
the  null  hypothesis  filter  path  in  Figure  2-1.  Of  course,  the 
detection  statistic  and  decision  rule  are  modified  appropriately. 

2 . 2  Scalar.  Two-Dimensional  Representation  of  the  Channel 
Output  Sequence 

The  channel  output  vector  sequence  {x(n)}  can  be  viewed  as  a 
scalar  2-D  sequence  (Dudgeon  and  Mersereau,  1984) .  Let  channel  J 
be  the  temporal  and  spatial  reference  for  the  array,  and  define 
the  following  association, 
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(2-1) 


Xj-k(n)  =  X(n>k) 


0  <  n  <  N-1;  0<k<J-1 


The  process  {x(n,k)}  has  a  scalar  2-D  power  spectrum  denoted  as 
{Sxx(fd,fs)},  and  a  scalar  2-D  auto-covariance  sequence  (ACS)  denoted 
as  {rxx(m,^)}/  where  (m,^)  is  the  lag  pair  corresponding  to  the 

frequency  pair  (fd,fs). 

The  association  defined  by  Equation  (2-1)  is  illustrated  in 
Figure  2-2;  this  association  was  adopted  herein  because  it  is 
consistent  with  the  relation  between  multichannel  1-D  and  scalar 

2-D  systems  demonstrated  by  Therrien  (1981) .  Alternatives  to 
Equation  (2-1)  can  be  defined,  such  as  {xk+1(n)  =  x(n,k)  I  n  =  0,  1 . N-1;k 

=  0,  1 ,  .  .  .  ,  J-1 } .  This  alternative  definition  corresponds  to  the 
MATLAB  software  default  array  definition  convention. 

A  block  diagram  equivalent  to  that  of  Figure  2-1  can  be 
defined  using  the  scalar  2-D  process  {x(n,k)}  as  its  input  (see 
Section  3.3).  In  such  a  diagram  the  functionality  of  the  blocks 
remains  the  same,  but  the  specific  content  of  each  block  must 
accommodate  the  differences  in  model  class  and  implementation 
aspects . 

As  stated  earlier,  the  time  series  subclass  of  2-D,  linear, 
shift-invariant  parametric  models  was  selected  to  represent  the 
channel  output  under  each  of  the  two  hypotheses.  The  work  in 
Phase  I  addressed  two  important  issues  beyond  the  work  of  Michels 
(1991) .  First,  the  channel  output  is  represented  as  a  2-D  scalar 
random  process,  with  the  attendant  increase  in  modeling  degrees- 
of-freedom  (independent  dynamic  and  static  modeling  capability 
along  each  axis) .  Second,  the  complete  time  series  model  class 
(MA,  AR,  and  ARMA  models)  was  considered,  with  emphasis  on  ARMA 
models  due  to  their  higher  degree  of  modeling  capacity. 
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2 . 3  Two-Dimensional  Time  Series  Models 

Model-based  detection  as  considered  herein  is  predicated  on 
the  representation  of  the  channel  output  process  (x(n,k)}  as  the 
output  of  a  2-D  time  series  model  driven  by  white  noise. 
Furthermore,  the  time  series  model  output  is  corrupted  by  additive 
white  noise.  To  be  precise,  such  a  representation  is  approximate 
in  the  case  of  radar  data.  Nevertheless,  in  practice  time  series 
models  have  been  shown  to  provide  a  good  fit  to  radar  data;  this 
point  is  underscored  for  2-D  models  by  the  results  in  Section  4.0. 

As  stated  previously,  the  2-D  process  {x(n,k)}  is  assumed  to  be 
represented  as 
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(2-2) 


x(n,k)  =  y(n,k)  +  w(n,k) 


where  y(n.k)  is  the  output  of  a  linear,  shift-invariant,  2-D  time 
series  model,  and  w(n,k)  is  a  2-D,  zero-mean,  Gaussian-distributed, 
white  noise  process.  The  processes  {y(n,k)}  and  {w(n,k)}  are 
independent . 

A  2-D  ARMA(P,Q)  process  {y(n,k)}  is  defined  as 
P  P  Q  Q 

(2-3)  y(n,k)  =  -  ]T  Xa(id>is)y(n-id,k-is)  +  2  2b(id,is)u(n-id,k-is) 

id=1  's=1  id=0  's=0 


where  the  input  {u(n,k)}  is  a  2-D,  zero-mean,  Gaussian-distributed, 
white  noise  process,  and  {a.(id,is)}  and  {b(id,is)}  are  complex-valued, 

constant  coefficients  referred  to  as  the  AR  and  the  MA  parameters, 
respectively.  When  P  =  0  Equation  (2-3)  is  an  MA(Q)  process, 
whereas  when  Q  =  0  Equation  (2-3)  is  an  AR(P)  process.  The  transfer 
function  associated  with  model  (2-3)  is  obtained  as  the  2-D  Z- 
transform  of  Equation  (2-3).  That  is,  the  transfer  function  is  a 
2-D,  complex- valued,  rational  function  of  the  form 


(2-4)  T(zd,zs) 


b(qd,qs)zdqdzsqs 
B(zd,zs)  ^  qd=oqs=o _ 

A(zd,zs)  v  v  /  \  -pd 

X  Xa(Pd-Ps)zd  z 


Pd  7-Ps 

S 


Pd=°  Ps=° 


(2-5)  a(0,0)  =  1 


with  complex-valued  variables  Zd  and  Zs,  and  where  A(zd,  zs)  and 
B(Zd-Zg)  are  the  AR  and  MA  scalar,  2-D  polynomials,  respectively. 
As  indicated  in  Equation  (2-5),  the  leading  coefficient  of  A(Zd,Zg) 

is  unity.  This  follows  from  Equation  (2-3),  and  is  the  2-D 
version  of  a  monic  polynomial  in  1-D.  And  the  2-D  frequency 
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response  of  model  (2-3),  which  is  denoted  as  T(fd,fs),  is  obtained  by- 
restricting  the  complex  variables  zd  and  Zs  to  the  unit  surface, 

(2-6)  zd  =  ej27lfd 

(2-7)  Zs  =  ej27lfs 

with  fd  and  fs  the  normalized  temporal  (Doppler)  and  spatial 
frequencies,  respectively. 

In  the  context  of  MIBDA  processing  for  surveillance  phased 
array  radar  systems,  modeling  involves  determination  of  the 
parameters  {a(id,is)}  and  {b(id,is)}  to  represent  a  moving  target,  jamming 

interference,  and  ground  clutter.  These  channel  output  components 
must  be  modeled  either  individually,  or  in  any  one  of  the  possible 
combinations,  as  dictated  by  the  scenario  conditions  and  the 
processor  architecture.  For  PAMF  processing  only  ground  clutter 
and/or  jamming  interference  must  be  modeled.  In  Phase  I,  emphasis 
was  placed  in  modeling  ground  clutter  since  it  drives  the 
requirements  for  a  2-D  model.  Having  identified  the  model 
parameters  which  represent  the  process,  the  associated  2-D 
whitening  (inverse)  filter  is  available  directly.  Specifically, 
given  the  channel  output  2-D  sequence,  {x(n,k)},  the  whitening  filter 
residual  sequence,  {v(n,k)},  is  obtained  as 

Q  Q  P  P 

(2-8)  v(n,k)  =  -  ^  ^b(id,is)v(n-id,k-is)  +  ]Ta(id,is)  x(n-id,k-is) 

's"1  *s=0 

If  model  (2-2) -(2-3)  represents  the  channel  output  exactly  and  the 
noise  {w(n,k)}  is  zero,  then  {v(n,k)}  is  a  true  innovations  sequence. 

For  any  other  conditions,  the  filter  residual  approximates  an 
innovations.  Notice  that  the  whitening  filter  is  also  an  ARMA 
system,  but  with  the  roles  of  the  AR  and  MA  coefficients  reversed. 
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A  2-D  model  of  the  form  in  Equation  (2-3)  is  referred  to  as  a 
first-quadrant  system  since  only  values  in  the  first  quadrant  of 
the  input  and  output  2-D  planes  (except  for  initial  conditions) 
are  used  to  generate  the  system  output.  Model  (2-3)  is  causal,  in 
loose  analogy  with  the  1-D  case,  since  only  past  outputs  and 
present  and  past  inputs  are  used  to  generate  the  present  output . 
In  the  phased  array  radar  space-time  problem  causality  along  the 
time  axis  is  an  inherent  feature  of  the  channel  output.  In 
contrast,  the  issue  of  causality  along  the  spatial  axis  is  not 
clear  because  all  channels  generate  an  output  at  each  temporal 
sampling  instant.  Model  (2-3)  is  also  recursively  computable 
(Dudgeon  and  Mersereau,  1984),  which  simplifies  hardware 
implementation.  Other  region  of  support  options,  such  as  the  non- 
symmetric  half  plane  (NSHP)  ,  are  of  interest  and  will  be 
considered  in  Phase  II. 

Two-dimensional  time  series  models  have  several  features 
quite  distinct  from  their  1-D  counterparts.  To  begin,  2-D  models 
offer  more  dynamic  as  well  as  static  modeling  degrees-of- freedom. 

On  the  negative  side,  most  2-D  polynomials  are  not  factorable. 
Thus,  polynomials  A(zd,zs)  and  B(zd,zs)  identified  by  an  algorithm  are 

likely  to  be  unf actorizable .  This  complicates  key  issues  such  as 
stability  determination.  In  2-D  systems,  poles  and  zeros  can 
occur  as  functions  rather  than  as  an  isolated  point  (Dudgeon  and 
Mersereau,  1984).  Also,  causality  and  region  of  support  issues 
often  present  a  variety  of  alternatives,  in  contrast  with  single 
options  in  1-D  systems.  All  cases  considered  in  Phase  I  resulted 
in  a  stable  2-D  model. 

Estimation  of  the  parameters  for  the  2-D  scalar  time  series 
models  has  been  addressed  by  several  authors  (see  [Marple,  1987] , 
[Dudgeon  and  Mersereau,  1984] ,  [Therrien,  1986]  and  the  references 
therein) .  For  AR  models  most  algorithms  require  generation  of  the 
ACS  of  the  process,  and  the  algorithms  are  extensions  of  the  1-D 


13 


case.  In  contrast  with  the  1-D  case,  for  2-D  systems  the  maximum 
entropy  method  (MEM)  is  not  equivalent  to  the  AR  method.  Further, 
the  MEM  parameters  are  obtained  via  optimization  procedures,  with 
their  attendant  convergence  and  other  such  difficulties.  Most  2-D 
model  identification  algorithms  are  formulated  in  a  stochastic 
framework  and  require  the  ACS  (in  practice,  only  an  estimate  of 
the  ACS  is  available)  .  An  exception  is  the  2D-LS-FD  method 
developed  recently  by  Mikhael  and  Yu  (1994).  The  2D-LS-FD 
algorithm  is  a  least-square  method,  as  its  name  implies,  and  is 
the  algorithm  selected  in  Phase  I. 

A  generalized  version  of  the  2-D  ARMA(P.Q)  process  can  be 
defined  by  allowing  distinct  values  for  the  model  order  along  each 
dimension.  Specifically,  a  2-D  ARMA(P,Q)  process  {y(n,k)}  is  defined 
as 

Pd  ps  Qd  Qs 

(2-9)  y(n,k)  =  -  £a(id,is)y(n-id,k-is)  +  ^b(id,is) u(n-id,k-is) 

id=1  is=1  id=0  is=0 

where  the  input  process  and  system  parameters  are  as  defined 
previously,  and  P  =  [Pd,Ps],  Q  =  [Qd,Qs].  The  transfer  function  and 

residual  for  model  (2-9)  follow  trivially.  A  system  of  the  form 
(2-9)  is  a  more  relevant  model  for  the  STAP  application  because 
the  data  duration  along  the  temporal  axis  is  usually  longer  than 
the  data  duration  along  the  spatial  axis  (the  data  matrix  is 
rectangular) . 
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3.0  TWO-DIMENSIONAL,  LEAST-SQUARE,  FREQUENCY -DOMAIN  (2D- 

LS-FD)  ALGORITHM  FOR  MIBDA  PROCESSING 

The  2D-LS-FD  method  is  based  on  a  modification  to  the  output 
error  model  method,  and  is  an  extension  of  prior  work  for  the  1-D 
case.  The  output  error  model  method  is  a  model  matching  approach, 
which  inherently  invokes  least-squares  optimization.  Such  a 
formulation  is  in  contrast  with  the  stochastic  formulation  common 
to  most  model  identification  algorithms.  The  method  can  be  viewed 
also  as  a  direct  approach,  since  the  algorithm  operates  on  the 
channel  output  data  directly,  without  the  requirement  to  estimate 
the  channel  output  ACS.  A  key  feature  of  the  2D-LS-FD  method  is 
that  the  parameters  of  the  2-D  polynomial  model  are  calculated  by 
solving  a  linear  set  of  equations  in  the  frequency  domain. 
Another  important  feature  of  the  2D-LS-FD  method  is  that  an 
indication  of  model  order  is  available  by  the  dimension  of  a 
subspace,  which  translates  into  verifying  the  rank  of  a  matrix. 
Additionally,  the  2D-LS-FD  method  admits  alternative  algorithmic 
variations,  which  provides  a  wealth  of  possible  approaches  to 
model  the  airborne  surveillance  phased  array  radar  return.  Due  to 
these  features,  the  2D-LS-FD  method  was  the  baseline  algorithm  for 
Phase  I . 

Originally,  the  2D-LS-FD  method  was  developed  for  real-valued 
data,  and  for  a  square  ARMA  model,  Equation  (2-3),  since  these  are 
the  conditions  that  apply  in  image  processing  applications  such  as 
data  compression  and  noise  cancellation.  In  Phase  I  the  method 
was  extended  to  handle  the  complex-valued  case  and  the  rectangular 
ARMA  model,  Equation  (2-9).  Both  extensions  turned  out  to  be 
straightforward.  The  main  issue  in  the  extension  to  complex¬ 
valued  data  is  to  insure  that  the  rules  of  complex-valued 
operations  (inner  product;  2-norm;  etc.)  are  satisfied.  A  MATLAB- 
based  software  implementation  of  a  fast  algorithm  was  developed 
also  in  Phase  I. 
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A  MATLAB  set  of  routines  was  generated  to  model  the  complex¬ 
valued  version  of  the  2D-LS-FD  algorithm.  These  routines  were 
tested  and  validated  using  several  test  cases,  including  an 
ARMA (1,1)  system,  and  a  signal  with  a  power  spectrum  consisting  of 
the  sum  of  three  Gaussian-shaped  peaks.  Once  the  software  was 
validated,  it  was  exercised  to  model  and  to  whiten  simulated 
airborne  surveillance  phased  array  radar  data  generated  with  the 
SSC-generated  simulation  for  a  variety  of  system  parameter  values 
and  scenario  conditions.  The  SSC  simulation  is  described  in  Vol. 
II  of  the  report  by  Roman  and  Davis  (1997)  . 

Model  order  determination  is  a  fundamental  issue  for  all 
model  identification  algorithms,  including  the  2D-LS-FD  method. 
The  approach  postulated  by  SSC  and  pursued  in  Phase  I  is  based  on 
the  characteristics  of  the  residual  generated  using  the  inverse 
(whitening)  system  of  the  identified  model. 

3.1  2D-LS-FD  Method  for  2-D  Model  Identification 


Consider  the  output  error  model  configuration  defined  in 
Figure  3-1,  which  is  presented  in  2-D  Z-transform  domain  (all 

functions  are  complex-valued).  In  Figure  3-1,  the  2-D  function 
Tu(zd,zs)  represents  an  unknown  system,  T(zd,zs)  represents  an  ARMA 

model  ( 2-9 ) ,  U(zd,zs)  is  the  2-D  driving  white  noise  function,  and 
E0(zd,zs)  is  the  2-D  error  function.  It  is  important  to  recognize 

that  the  unknown  system  is  unrestricted  (that  is,  it  can  be 
different  from  an  ARMA) .  In  equation  form,  the  output  error  is 

(3-1)  E0(zd,zs)  =  U(zd,zs)  [Tu(zd,zs)  -  T(zd,zs)] 

and  the  2-norm  of  the  output  error  restricted  to  the  unit  surface 
(Equations  (2-6)  and  (2-7))  is  defined  as 


16 


(3  2)  E0 (f^ , fs )|  2 


0.5  0.5 


1 1/2 


J  J|E0(fd,f8)|2dfddfs 
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Minimization  of  the  2-norm  in  Equation  (3-2)  with  respect  to  the 
ARMA  model  parameters  generates  a  model  which  is  an  optimum 

approximation  to  the  unknown  system  for  the  selected  model  order. 
The  function  U(fd,fs)  serves  as  a  weighting  term;  as  such,  Equation 

(3-2)  represents  a  weighted  2-norm  minimization. 

Most  methods  available  for  minimizing  (3-2)  are  nonlinear, 
which  presents  computational  problems  and  convergence  issues. 
Recently,  however,  Mikhael  and  Yu  (1994)  have  developed  a 
formulation  which  leads  to  a  linear  method.  Their  formulation  is 
based  on  the  variation  of  Figure  3-1  shown  in  Figure  3-2  for  the 
case  of  an  ARMA  model,  and  is  referred  to  as  the  equation  error 
model  method.  In  this  report  the  equation  error  is  denoted  with 
subscript  "m"  in  order  to  emphasize  the  difference  between  the  two 
formulations.  Notice  that  all  2-D  functions  in  Figure  3-2  are 
restricted  to  the  unit  surface  in  the  frequency  domain.  In 
equation  form,  the  equation  error  of  Figure  3-2  is 

0-3 )  Em(fd,y  =  u(fd,y  [T0(fd,y  A(fd,y  -  B(fdy] 

Minimization  of  the  2-norm  of  the  equation  error  Em(fd,fs)  leads  to  a 
linear  set  of  equations  in  the  unknown  parameters  {a(id,is)}  and 
{b(id,is)}.  The  linear  set  of  equations  has  a  solution  provided  the 
number  of  equations  is  more  than  or  equal  to  the  number  of  unknown 

parameters,  and  provided  that  the  parameter  set  is  constrained 
appropriately.  In  this  modified  formulation  the  function  U(fd,fs) 

still  plays  the  role  of  a  weighting  term. 
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Figure  3-1.  Output  error  model  formulation  for  2-D  system 

modeling. 


u(fd,fs) 


Figure  3-. 


W.) 


Equation  error  model  formulation  for  input-driven  2-D 
system  modeling. 


3 . 2  Space-Time  Data  Modeling  Using  the  2D-LS-FD  Method 


In  both  formulations  discussed  above  (Figures  3-1  and  3-2) 
the  unknown  system  is  available  to  be  driven  by  the  input  function 
{u(n,k)>/  and  the  input  function  is  known.  Neither  of  these  two 
conditions  is  true  for  modeling  the  space-time  phased  array  radar 
signal  of  interest  herein.  Specifically,  in  the  STAP  channel 
output  modeling  case  only  the  output  of  the  "unknown  system"  is 
available.  Here  the  quotation  marks  emphasize  the  fact  that  the 
channel  output  is  not  the  output  of  a  particular  ARMA  system,  but 
rather,  that  it  is  to  be  approximated  as  such.  In  the  model 
matching  sense,  the  channel  output  sequence  is  the  output  of  an 
unknown  model  with  input  {u(n,k)}.  If  the  input  is  the  unit  sample 
sequence,  u(n,k)  =  5(n,k) ,  then  {x(n,k)}  can  be  viewed  as  the  impulse 

response  of  the  unknown  system, 

(3-4 )  x(fd,fs)  =  Tu(fd,y  U(fd,y  =  Tu(fd,y 

In  Equation  (3-4)  the  second  equality  is  valid  because  the 
discrete  Fourier  transform  (DFT)  of  the  unit  sample  sequence  is 

(3-5)  U(fd,fs)  =  1  -0.5  <  fd  <  0.5;  -0.5<fs<0.5 

For  this  input  condition,  the  equation  error  is  expressed  as 

0-6)  E(fd,y  =  x(f„,fs)  A(fd,y  -  B(fd,y 

This  equation  error  formulation  is  denoted  without  subscript  in 
order  to  emphasize  the  difference  with  the  previous  cases.  Notice 

the  similarity  in  form  between  the  expressions  for  the  equation 
errors  Em(fd,fs)  and  E(fd,fs).  A  block  diagram  for  the  equation  error 

formulation  (3-6)  is  illustrated  in  Figure  3-3,  where  the  dashed 
lines  represent  the  condition  that  the  unknown  system  is 
unavailable  and  thus  cannot  be  driven  by  the  input  sequence. 
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U(fd,fs) 


E(fd,fs) 


Figure  3-3.  Equation  error  model  formulation  for  modeling  the 
space- time  phased  array  radar  process. 


Minimization  of  the  2-norm  of  the  equation  error  E(fd,fs)  with 

respect  to  the  unknown  ARMA  coefficients  is  a  finite-dimensional, 
linear-space  minimization  problem.  Mikhael  and  Yu  (1994)  applied 
the  Orthogonality  Theorem  to  the  problem,  which  leads  to  a  linear 
set  of  equations  in  the  unknown  parameters  {a(id,is)}  and  {b(id,is)}. 

Such  equations  have  a  solution  provided  the  number  of  equations 
exceeds  the  number  of  unknown  parameters.  Furthermore,  it  is 

necessary  to  place  appropriate  constraints  in  order  to  avoid  the 
trivial  solution,  A(fd.f  s)  =  0  and  B(fd,fs)  =  0 .  One  such  constraint  is  to 

set 

(3-7a)  A(fd,fs)  =  1  +  A0(fd,fs)  =  1+£ 

Pd=°  Ps=° 


(3-7b)  a(0,0)  =  0 

This  constraint  is  equivalent  to  the  definition  of  A(fd,fs)  given  in 
Equations  (2-4)  and  (2-5)  .  Also,  notice  that  the  monic  feature  of 
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the  2-D  polynomial  A(fd,fs)  is  retained.  Substitution  of  this 
constraint  into  Equation  (3-6)  results  in 

( 3  -  8 )  E(fd,fs)  =  X(fd,fs)  -  [B(fd,fs)  -  X(fd,fs)  A0(fd,fs)]  =  X(fd,fs)  -  G(fd,fs) 

with  the  2-D  function  G(fd,fs)  defined  implicitly  as 
(3-9)  G(fd,fs)  =  B(fd,fs)  -  X(fd,fs)  A0(fd,fs) 

Function  G(fd,fs)  is  an  auxiliary  function  introduced  to  simplify  the 

analysis.  Minimizing  the  2-norm  of  the  equation  error  with 

respect  to  the  unknown  coefficients  corresponds  to  minimizing  the 
difference  X(fd,fs)  -  G(fd,fs) .  Computational  advantages  accrue  because 

the  problem  is  formulated  in  the  frequency  domain,  and  the 
resulting  solution  is  optimal  in  the  least-square  sense  (Mikhael 
and  Yu,  1994)  .  Once  the  system  model  parameters  have  been 

identified,  the  inverse  (whitening)  system  is  obtained  as  the 
inverse  of  Equation  (2-4)  ,  since  T(fd,fs)  is  a  scalar  function. 

This  method  is  related  to  that  proposed  by  Shanks  et  al . 
(1972)  for  the  design  of  2-D  infinite  impulse  response  filters. 
Shanks'  method,  however,  is  formulated  in  the  space-time  domain 
(as  opposed  to  the  frequency  domain) ,  and  its  algorithmic 
implementation  is  different. 

3 . 3  Two-Dimensional  MIBDA 


In  the  approach  pursued  in  Phase  I,  the  MIBDA  in  Figure  2-1 
is  applicable  for  2-D  processing  also,  with  the  modification  that 
a  2-D  ARMA  whitening  filter  is  generated  for  each  of  the  two 
hypotheses  in  the  MIBDA.  However,  2-D  convolution  in  the  time 
domain  is  a  computationally-intensive  operation,  so  an 
architecture  with  frequency-domain  convolution  is  preferred.  Such 
an  architecture  is  presented  in  Figure  3-4  for  the  case  of  off- 
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line  filter  design.  A  similar  architecture  can  be  defined  for  the 
case  of  on-line  filter  design.  In  the  architecture  of  Figure  3-4, 
the  Pre-Processor  block  can  include  a  variety  of  operations. 
Several  such  functions  were  identified  in  Phase  I,  and  evaluated 
to  different  degrees  of  detail.  The  most  important  operations  for 
the  Pre-Processor  block  are  presented  next,  in  the  preferred  order 
of  application  to  the  data. 

•  Spatial  nulling  of  the  1-D  vector  sequence  (this  eliminates 
the  effects  of  barrage-type,  point-location  jammers,  which 
simplifies  the  data  structure  in  the  frequency  domain) . 

•  1-D  vector  to  2-D  scalar  transformation  (this  is  only  a 
conceptual  operation  since  the  storage  locations  of  the 
data  remain  unchanged) . 

•  Application  of  a  data  window  (this  reduces  the  effects  of 
sidelobe  leakage,  at  the  expense  of  a  loss  in  resolution, 
which  results  in  broader  spectral  features) . 

In  a  given  configuration,  only  a  subset  of  these  operations  may  be 
applied  to  the  data.  Optimization  of  the  pre-processing  functions 
will  be  carried  out  in  the  proposed  Phase  II.  This  includes 
analyses  such  as  performance  variation  as  a  function  of  window 
type  applied. 


te(n)} 


Figure  3-4.  Multichannel  innovations -based  detector  architecture 

for  2-D  data  representation. 
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The  2-D  FFT  operation  can  be  followed  by  an  averaging  step, 
which  is  not  shown  in  Figure  3-4.  Averaging  the  spectral 
components  over  several  realizations  generates  a  more  accurate 
estimate  of  the  frequency-domain  representation.  This  statement 
is  applicable  during  hypothesis  filter  design,  as  well  as  in 
processor  architecture  design. 

It  is  important  to  note  that,  with  respect  to  pre-processing 
and  spectral  averaging,  the  MIBDA  design  must  be  compatible  with 
the  procedure  utilized  during  hypothesis  filter  design.  For 
example,  if  a  data  window  of  a  specific  type  is  used  in  the  design 
of  the  hypothesis  filters,  then  the  same  data  window  must  be 
applied  in  the  MIBDA  implementation.  Likewise,  if  the  hypothesis 
filters  are  designed  with  spectral-domain  averaging  over  K 
realizations,  then  the  MIBDA  must  include  spectral-domain 
averaging  over  K  realizations  also. 

Consider  now  the  2-D  inverse  FFT  block  in  Figure  3-4,  which 
includes  implicitly  the  2-D  scalar  to  1-D  vector  transformation. 
This  block  is  necessary  if  the  hypothesis  filter  design  criteria 
and/or  the  detection  criteria  and  decision  rules  are  to  be  applied 
in  the  space-time  domain,  as  opposed  to  the  frequency  domain.  In 
a  prior  program  (Roman  and  Davis,  1997)  SSC  developed  a  set  of 
hypothesis  filter  design  criteria  that  are  applied  in  the  space- 
time  domain  for  the  1-D  vector  MIBDA.  Those  criteria  have  been 
extended  in  Phase  I  to  the  2-D  scalar  case,  but  remain  to  be 
implemented  in  software.  Similarly,  the  likelihood  ratio  detector 
used  in  the  1-D  vector  MIBDA  can  be  extended  to  the  2-D  scalar 
case  and  applied  in  the  space- time  domain.  Thus,  based  on  the 
criteria  considered  thus  far,  the  inverse  FFT  step  is  required. 
However,  the  alternative  of  frequency-domain  detection  is 
computationally  attractive,  and  will  be  investigated  in  the 
proposed  Phase  II. 


3 . 4  Model  Order  Determination 


Determination  of  model  order  is  a  decision  required  of  all 
model  identification  algorithms  in  applications  where  the  true 
order  of  the  system  generating  the  channel  output  data  is  unknown, 
or  where  the  process  to  be  modeled  is  generated  by  a  system  which 
belongs  to  a  different  model  class.  In  the  second  case  the  model 
obtained  is  a  "representation  model,"  as  opposed  to  a  "physical 
model"  (a  model  based  on  accurate  analyses  of  the  underlying 
physical  processes) .  Model  order  determination  is  always  a 
difficult  problem,  and  the  approaches  are  often  heuristic,  at 
least  in  part.  The  approach  selected  herein  is  based  on  the 
output  power  (variance)  of  the  model  inverse  (whitening  filter) , 
since  the  best-fitting  model  would  generate  the  most  whitening  (in 
the  sense  of  whitening  the  ensemble,  not  just  a  particular 
realization) .  Such  a  criterion  has  to  be  applied  judiciously, 
however,  because  it  can  lead  to  over-parameterized  models  since 
the  model  with  the  most  possible  degrees-of-f reedom  (number  of 
parameters  to  be  identified)  is  likely  to  whiten  a  specific 
realization  the  most.  Thus,  for  a  criterion  based  on  the  residual 
power,  the  model  order  sought  is  that  value  at  which  the  plot  of 
residual  power  versus  model  order  indicates  diminishing  returns 
for  further  increases  in  model  order  (such  a  value  is  indicated  by 
a  "knee"  in  the  curve) . 

For  a  zero-mean  2-D  residual  sequence  of  the  form  in  Equation 
(2-8) ,  the  power  (variance)  is  estimated  as 


(3-10) 


6?  = 


1 

JN 


N-1  J-1 

T.  y>(n.K)f 


n=0  k=0 


In  the  context  of  model  order  determination  adopted  herein,  the 
power  estimate  Gy  is  plotted  as  a  function  of  the  order  of  the 
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model  used  to  generate  the  residual  sequence,  and  the  model  order 
is  selected  as  the  value  at  which  the  reduction  in  residual  power 
is  small  (in  relation  to  prior  reductions)  for  a  unit  increase  in 
model  order. 

A  better  criterion  is  one  which  considers  also  the  effect  of 
the  given  number  of  2-D  data  points  (JN)  and  the  number  of 
parameters  being  estimated  for  the  model.  Akaike  (1969;  1970) 

proposed  a  simple  criterion  for  AR  model  order  determination 
referred  to  as  the  final  prediction  error  (FPE)  which  accounts  for 
these  two  effects.  In  the  FPE  criterion  the  residual  power 
estimate  obtained  using  the  Yule-Walker  AR  parameter  estimation 
algorithm  is  modified  by  a  multiplicative  factor  which  is  a  non¬ 
decreasing  function  of  the  number  of  data  points  and  the  number  of 

parameters  being  estimated  for  the  model  order  being  evaluated. 
For  a  2-D  ARMA(Pd,Ps,Qd,Qs)  model  and  a  JxN  data  matrix,  the 

multiplicative  factor  corresponding  to  the  one  used  by  Akaike  is 
of  the  form 

JN  +  M 

(3-11)  f(M)  =  — — — — —  0  <  M  <  JN 

JN-M 

where  M  is  the  number  of  parameters  being  estimated, 

(3-12)  M  =  (Pd  +  1)(PS  +  1)  +  (Qd  +  1)(QS  +  1)  -  1 
And  for  the  special  case  where  Qd  =  Pd  and  QS  =  PS, 

(3-13)  M  =  2  (Pd  +  1)(PS  +  1)  -  1 

Notice  that  the  function  f(M)  is  non-decreasing.  Thus,  Akaike ' s 
FPE  criterion  takes  into  consideration  the  following  two 
conditions:  (a)  as  the  model  order  increases,  the  fit  of  the  model 
to  the  given  data  improves  (the  residual  power  should  decrease 
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monotonically  with  increasing  model  order) ;  and  (b)  as  the  model 
order  increases,  the  fit  of  the  model  to  an  independent 
realization  of  the  same  stochastic  process  becomes  worse  on  the 
average  (modeling  degrees-of-f reedom  beyond  those  necessary  to 
model  a  process  are  allocated  by  an  identification  algorithm  to 
modeling  features  present  only  in  the  realizations  used  to 
identify  the  model) . 

Akaike's  FPE  criterion  has  been  the  subject  of  multiple 
studies,  many  of  which  resulted  in  a  modification  to  the  criterion 
expression.  In  most  cases  the  FPE  has  provided  a  good  estimate  of 
the  true  model  order.  The  extension  of  the  FPE  to  the  2-D  case 
will  be  evaluated  in  the  proposed  Phase  II. 


4 . 0  AIRBORNE  SURVEILLANCE  PHASED  ARRAY  RADAR  APPLICATION 


Airborne  surveillance  for  moving  target  detection  using 
phased  array  radar  systems  is  one  of  the  major  technology  thrusts 
at  RL.  Consequently,  the  emphasis  devoted  to  that  application  in 
Phase  I  reflects  that  significance.  The  MATLAB  software 
implementation  of  the  2D-LS-FD  algorithm  was  exercised  with 
simulated  multichannel  data  generated  using  the  SSC-developed 
software  simulation  package  (Roman  and  Davis,  1997),  and  selected 
results  are  presented  herein. 

The  SSC  airborne  surveillance  phased  array  radar  software 
model  (Vol.  II  of  [Roman  and  Davis,  1997])  is  defined  for  a  side¬ 
looking,  linear',  phased  array  radar  configuration,  with  azimuth 
scanning  capability  of  ±90  deg  from  boresight  (the  effect  of 
multiple  elevation  channels  beamformed  into  a  single  channel  is 
included  in  the  model) .  Each  component  (target;  noise;  jammer; 
clutter)  is  modeled  independently  of  the  others.  The  structure 
selected  to  model  each  component  allows  generation  of  the  true 
covariance  matrix  sequence,  as  well  as  generation  of  independent 
statistical  data  realizations.  Moving  targets  are  modeled  as 
point  targets  in  azimuth  and  elevation,  and  each  moving  target  is 
modeled  as  the  output  of  a  first-order  linear  system  in  state- 
space  form.  This  allows  for  a  wide  range  of  conditions  (including 
Swerling  Cases  1  through  4  as  well  as  the  deterministic  case) ,  in 
a  simple  and  general  format.  With  this  structure  Gs  targets  are 
grouped  as  a  Gsth-order  linear  system.  Barrage- type  jammers  are 
modeled  as  point  sources  of  Gaussian-distributed  broadband  noise 
in  azimuth  and  elevation.  Receiver  noise  is  modeled  as 
uncorrelated  in  time  (pulse-to-pulse)  as  well  as  in  space 
( channel -to-channel) .  Ground  clutter  is  modeled  as  a  large  number 
of  statistically- independent  clutter  patches,  with  each  clutter 
patch  consisting  of  a  large  number  of  individual  re-radiators. 


Table  4-1  lists  the  baseline  simulation  parameter  values  used 
to  obtain  most  of  the  modeling  and  filtering  (whitening)  results 
presented  in  this  section.  Clutter-to-noise  ratio  (CNR)  for  both 
sets  of  the  baseline  conditions  (crab  angle  values  of  y=0  deg  and 
y  =  20  deg)  is  47.75  dB.  Notice  that  jammers  are  not  included  in 
these  results;  this  is  because  the  emphasis  of  the  simulation- 
based  analyses  has  been  for  ground  clutter  in  receiver  noise  since 
it  presents  the  most  significant  modeling  challenge.  Also,  in  one 
of  the  alternative  configurations  jammers  are  canceled  spatially 
in  the  pre-processing  step.  Modeling  and  whitening  runs  with 
jammers  present  have  been  made,  and  the  results  are  similar  to 
those  for  clutter  only.  The  model  structure  analysis  results 
shown  herein  were  obtained  for  parameter  values  and  simulation 
conditions  different  from  those  in  the  table;  the  actual  parameter 
values  used  in  such  cases  are  stated  clearly.  Several  cases  were 
run  also  with  the  factored  form  of  the  classical  optimal  joint- 
domain  method  (Brennan  and  Reed,  1973),  in  order  to  compare  the 
whitening  capability  of  the  optimum  method  with  that  of  the  2D-LS- 
FD  method.  Those  results  show  that  the  2D-LS-FD  method  provides 
significantly  more  whitening  capability  than  the  classical  optimal 
joint-domain  method  with  a  limited  support  size  for  the  sample 
covariance  matrix. 

Alternative  modeling  configurations  of  the  2D-LS-FD  method 
(ARMA;  AR;  MA)  were  evaluated  in  the  context  of  the  airborne 
surveillance  radar  array  STAP  application,  with  ARMA  modeling 
providing  the  best  model  fit  and  whitening  performance.  This 
result  is  as  expected  since  the  ARMA(P,Q)  model  with  P  >  Q  provides 
a  larger  number  of  modeling  degrees-of-f reedom  than  an  AR(P)  or  an 
MA(P) ,  with  the  same  dynamic  realization  requirements. 
Furthermore,  the  ARMA(P,Q)  model  with  P>Q  offers  more  versatility 
than  an  AR(P+Q)  or  an  MA(P+Q) .  Simulation  results  also  indicate 
that  an  ARMA(P.P)  model  in  general  provides  better  modeling  than  an 
ARMA(R,Q)  with  R  <  Q  and  R  +  Q  <  2P . 
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PARAMETER 

TYPE 

PARAMETER  (UNITS) 

RADAR 

Number  of  linear  array  elements,  J 

Number  of  points  in  one  CPI,  N 

Number  of  elevation  axis  elements,  Je 

Array  mainbeam  azimuth  angle,  phiO  (deg) 

Peak  transmitted  power,  Pt  (kW) 

Pulse  (uncompressed)  duration,  Tu  (jisec) 

Pulse  repetition  frequency  (Hz)  j 

ARRAY 

Radiation  frequency,  fC  (MHz) 

SYSTEM 

Receiver  bandwidth,  fB  (MHz) 

Transmit  pattern  gain,  Go  (dB) 

Receive  element  gain,  Ge  (dB) 

Receive  element  backlobe  pattern  attenuation,  Gb  (dB) 

Noise  figure,  Fn  (dB) 

System  losses,  Ls  (dB) 

Transmit  pattern  array  option,  patopt 

SURVEILLANCE 

Platform  altitude,  Hp  (km)  j 

Platform  velocity,  Vp  (m/sec) 

SCENARIO 

Range  to  principal  ground  clutter  ring,  rc  (km) 

Aircraft  platform  crab  angle,  gamma  (deg) 

TARGET 

Narrowband  process  amplitude,  a 

Target  radial  velocity,  Vt  (m/sec) 

Target  azimuth  angle,  phit  (deg) 

Target  elevation  angle,  thetat  (deg) 

Signal-to-noise  ratio,  SNR  (dB) 

INTERFERENCE 

Jammer  azimuth  angle,  phii  (deg) 

Jammer  elevation  angle,  thetai  (deg) 

Jammer  power,  vari 

GROUND  CLUTTER 

Number  of  ground  patches  illuminated  by  mainbeam,  Nc 

ARRAY  NOISE 

Receiver  noise  power  per  channel,  varn 

SIMULATION 

Number  of  realizations  used  in  filter  design,  Nrd 

PARAMETERS 

Number  of  realizations  used  in  filter  evaluation,  Nre 

VALUE 


UNIFORM 


9 


5 


130 


0;  20 
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4 . 1  Conventional  Optimal  Joint -Domain  Method 

Consider  the  JN-element  data  vector  Xk  formed  by  concatenating 
the  N  J-element  vectors  {x(n)}  of  the  kth  range  cell  in  a  coherent 
processing  interval  (CPI) ,  and  let  denote  the  JNxJN  sample 
covariance  matrix  estimate  obtained  by  averaging  over  K 
independent  realizations  of  the  data  vector  (K  range  cells,  not 
necessarily  contiguous) .  That  is, 


(4-1) 


x  xH 
-k-k 


As  is  well-known,  the  optimal  weight  vector  is  determined  as 
(4-2)  W0  =  ^S  =  (^'2)(r'2s) 


where  S  is  the  known,  JN-element  steering  vector.  The  optimal 
weight  vector  is  applied  to  the  data  in  order  to  generate  the 
detection  statistic,  r|,  as 

( 4-3  )  T|  =  wj  X  =  (§H  %iV2  ) (^'1/2  ) X 

Equations  (4-2)  and  (4-3)  exhibit  the  factored  form  of  the  optimal 

~  -1/2 

joint-domain  weight  vector:  the  %  operator  is  a  whitening 
filter  for  X,  and  the  SM  ^  operator  is  a  matched  filter.  This 
factored  form  is  illustrated  as  a  block  diagram  in  Figure  4-1. 
The  intermediate  variable  in  the  factored  form, 

(4-4)  Y0  =  X 

is  defined  herein  as  a  JN  -element  residual  vector.  This  residual 

vector  can  be  transformed  into  a  data  array  and  subsequently  into 
a  2-D  scalar  sequence,  {v0(n,k)}(  by  an  association  analogous  to  the 
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one  established  in  Figure  (2-1) .  Then,  the  residual  optimal 
conventional  2-D  scalar  residual  sequence  can  be  compared  for 
whiteness  with  the  2-D  scalar  residual  sequence  {v(n,k)}  generated  by 

the  2-D  whitening  filter  designed  based  on  the  2D-LS-FD  algorithm. 


Range  Cell 
Data  Vector 


V 


Residual 


WHITENING 

FILTER 

-ok  ^ 

MATCHED 

FILTER 

xV2 

sH  hi"2 

Decision 

Statistic 


Figure  4-1.  Conventional  optimal  joint-domain  method  in  factored 

form. 


The  level  of  whitening  can  be  assessed  using  the  whitening 
ratio  performance  measure,  which  is  defined  herein  as 


(4-5) 


WR  =  10  log10 


N-1  J-1 


X  Xiv(n'k)i: 


n=0  k=0 


-10  log 


10 


N-1  J-1 


X  Xix(n>k>f 


n=0  k=0 


This  measure  can  be  applied  also  to  a  conventional  optimal  joint- 
domain  residual  by  substituting  v0  (n.k)  in  place  of  v0(n,k) . 

4 . 2  Simulated  Ground  Clutter  Modeling  and  Filtering 

In  the  context  of  time  series  modeling,  selection  of  the 
model  structure  (AR  vs.  MA  vs.  ARMA)  must  be  carried  out  prior  to 
model  order  determination.  However,  the  two  issues  are  inter¬ 
related  in  a  natural  manner,  and  one  approach  to  analyze  model 
structure  for  a  given  application  is  to  evaluate  the  three  model 
types  as  a  function  of  model  order.  Such  an  analysis  was  carried 
out  in  Phase  I  for  the  baseline  conditions  of  Table  4-1,  with  the 
following  alterations:  J  =  N  =  16;  and  Nrcj  =  64  (number  of 
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realizations  used  in  filter  design,  which  becomes  the  number  of 

spectra  averaged  for  model  identification  with  the  2D-LS-FD 
algorithm).  Also,  Pd  =  Ps  =  P  and  Qd  =  Qs  =  Q  were  assumed.  Figure  4- 

2  shows  a  plot  of  the  2-D  residual  power  as  a  function  of  model 
order  for  the  y=0  deg  case  and  the  following  model  structures  fit 

to  the  ground  clutter  return:  (a)  an  MA(Q)  model  (dotted  line) ; 
(b)  and  AR(P)  model  (dash-dot  line) ;  (c)  and  ARMA(Q-3,Q)  model 

(dashed  line)  ;  and  (d)  and  ARMA(P,P)  model  (solid  line) .  Notice 
that  the  ARMA(P,P)  model  generates  the  smallest  residuals 
(considering  a  "smoothed"  fit  to  the  curves  shown) .  This  result 
is  appealing  in  the  context  of  MIBDA  for  the  application  of 
interest  because  the  whitening  filter  associated  with  an  ARMA  is 
also  an  ARMA. 

The  results  in  Figure  4-2  are  based  on  only  one  generated 
residual  sequence  for  each  plotted  point,  but  similar  power  values 
have  been  obtained  with  other  independent  runs.  Also,  the  same 
conclusion  has  been  obtained  with  other  sets  of  parameter  values. 
Therefore,  the  ARMA(P,P)  model  structure  was  selected  for  other 
analyses,  but  additional  model  structure  analyses  will  be  made  in 
the  proposed  Phase  II . 

Consider  now  the  ARMA(P.P)  model  structure  as  a  function  of 

model  order  for  the  baseline  conditions  in  Table  4-1,  including 
the  two  crab  angle  value  options  (y=0  deg  and  y=20  deg)  .  Two- 

dimensional  residual  power  as  a  function  of  model  order  for  an 
ARMA(P,P)  model  fit  to  the  ground  clutter  return  is  presented  in 
Figure  4-3  for  both  crab  angle  conditions.  Each  point  of  each 
curve  in  Figure  4-3  is  an  average  of  twenty  (20)  independent  runs, 
which  accounts  for  the  smoother  character  of  these  two  curves  in 
relation  to  the  curves  in  Figure  4-2 .  In  Figure  4-3  notice  that 
the  minimum  residual  power  is  obtained  with  P  =  6  for  both 
conditions.  Residual  power  with  P  =  6  is  25.2  dB  for  y=0  deg, 

corresponding  to  a  -40.0  dB  whitening  ratio,  and  residual  power 
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with  P  =  6  is  24.2  dB  for  y=20  deg,  corresponding  to  a  -41.0  dB 
whitening  ratio.  This  whitening  performance  is  close  to  the 
maximum  value  attainable,  the  bound  set  by  the  analytical  CNR  of 
47.75  dB. 

The  parameter  values  in  Table  4-1  were  used  to  generate  two 
sets  of  whitening  performance  results,  indexed  by  the  two  baseline 
values  of  crab  angle,  y  =  0  deg  and  y=20  deg.  Whitening  results 

(same  conditions  but  different  data  realizations)  were  obtained 
also  for  the  classical  optimal  joint-domain  method  in  factored 
form.  Those  results  are  discussed  next. 

A  3-D  plot  of  the  channel  output  power  spectrum  is  presented 
in  Figure  4-4  for  the  y=0  deg  condition.  This  plot  is  an  average 

of  20  independent  2-D  modified  periodograms  ("modified"  denotes 
that  a  data  window  is  applied  prior  to  the  FFT) .  For  comparison, 
the  power  spectrum  of  the  2-D  ARMA(6,6)  model  identified  from  the 
data  is  presented  in  Figure  4-5.  Notice  the  similarity  between 
these  two  spectral-domain  figures. 

Figure  4-6  presents  a  plot  of  the  real  part  of  the  1-D 
circular  ACS  of  the  channel  3  sequence  prior  to  filtering  for  y=0 

deg.  This  plot  is  in  sharp  contrast  with  Figure  4-7,  which  is  a 
plot  of  the  real  part  of  the  1-D  circular  ACS  of  the  channel  3 
residual.  The  filtered  sequence  is  indeed  approximately  white. 
Analogously,  Figure  4-8  presents  a  plot  of  the  real  part  of  the  1- 
D  circular  ACS  of  the  received  signal  versus  spatial  lag  (k)  across 
the  J  =  8  channels  at  the  time  instant  corresponding  to  n  =  16.  And 
Figure  4-9  presents  a  plot  of  the  1-D  circular  ACS  of  the  real 
part  of  the  residual  across  the  J  =  8  channels  at  the  time  instant 
corresponding  to  n=16.  Notice  that  the  unfiltered  sequence  is 
highly  correlated,  whereas  the  residual  is  significantly  less 
correlated  (it  is  difficult  to  assess  the  level  of  whiteness  for  a 
short-duration  sequence).  These  four  figures  correspond  to  y=0 
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deg  conditions.  Also,  similar  results  are  true  for  the  imaginary- 
part  of  the  respective  ACS. 

Figures  4-10  and  4-11  present  the  real  part  of  the  circular 
ACS  of  the  residual  as  generated  by  the  whitening  step  of  the 
optimal  joint-domain  method  (see  Equation  (4-4))  for  y=0  deg 

conditions.  These  two  figures  can  be  compared  with  Figures  4-7 
and  4-9,  respectively.  For  these  two  figures  the  sample 
covariance  matrix  was  generated  using  Nre  =  JN  =  512  independent  CPI 
realizations,  which  is  the  minimum  number  of  realizations  required 
to  obtain  a  full-rank  covariance.  This  number  is  25  times  larger 
than  the  number  of  independent  realizations  (Nrd  =  Nre  =  20)  averaged 
to  generate  the  2-D  frequency-domain  representation  used  to  design 
and  to  evaluate  the  ARMA(6,6)  model.  Clearly,  the  ARMA-based 
whitening  filter  yields  residuals  with  significantly  less 
correlation  than  the  residuals  of  the  optimum  joint-domain  method. 

The  whitening  ratios  calculated  for  the  conventional  optimal 
joint-domain  residual  are  -30.15  dB  and  -21.96  dB  for  the  y=0  deg 
and  y =  20  deg  conditions,  respectively.  These  figures  are 
significantly  less  than  those  for  the  2-D  ARMA(6,6)  residuals 
(which  are  -40.0  dB  and  -41.0  dB,  respectively).  This  discrepancy 
is  due  to  the  less-than-optimal  (although  large)  number  of 
independent  realizations  used  to  generate  the  JNxJN  sample 
covariance  matrix  estimate. 

Figures  4-12  through  4-19  present  the  same  type  of  plots  as 
in  Figures  4-4  through  4-11,  with  the  following  variations.  First 
and  foremost,  Figures  4-12  through  4-19  were  generated  for  crab 
angle  y=20  deg  conditions.  Second,  channel  4  is  selected  for  the 

temporal  lag  plots,  instead  of  channel  3.  Third  and  last,  the 
spatial  lag  plots  are  for  pulse  n  =  44,  instead  of  n  =  16.  Careful 
examination  of  the  results  presented  for  y=20  deg  conditions  re¬ 
enforces  the  observations  derived  from  the  earlier  figures. 
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model  order,  P  or  Q 


Figure  4-2.  Residual  power  as  a  function  of  model  order  for  2-D 
MA,  AR,  and  ARMA  models  (J  =  N=  16;  Nrcj  =  64;  y=0  deg). 


model  order,  P 

Figure  4-3.  Residual  power  as  a  function  of  model  order  for  a  2-D 
ARMA(P,P)  model  for  both  crab  angle  conditions . 
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Figure  4-4.  Channel  output  log  power  spectrum  for  y=0  deg. 
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Figure  4-5.  Two-D  ARMA(6,6)  log  power  spectrum  for  y=0  deg. 
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temporal  lag,  m 


Figure 


4-10.  Real  part  of  channel  3  joint-domain  residual 
circular  ACS  versus  temporal  lag  (y=0  deg). 


REAL  PART  OF  PULSE  16  CLASSICAL  RESIDUAL  CIRCULAR  ACS 


Figure  4-11.  Real  part  of  pulse  16  joint-domain  residual  circular 

ACS  versus  spatial  lag  (Y=0  deg). 
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Figure  4-12.  Channel  output  log  power  spectrum  for  y=20  deg. 
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Figure  4-13.  Two-D  ARMA(6,6)  log  power  spectrum  for  y=20  deg. 
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5 . 0  CONCLUSIONS  AND  RECOMMENDATIONS 


In  Phase  I  the  STAP  problem  for  airborne  surveillance  using 
phased  array  radar  systems  was  formulated  as  a  2-D  model-based 
detection  problem,  and  a  2-D  innovations-based  methodology  was 
formulated.  Additionally,  a  model  identification  configuration 
based  on  the  2D-LS-FD  method  was  developed,  and  the  capability  to 
model  the  received  channel  process  (with  emphasis  on  the  ground 
clutter  process)  with  a  2-D  time  series  models  was  established. 
The  2D-LS-FD  method  (Mikhael  and  Yu,  1994)  was  selected  for  this 
study  because:  (a)  model  coefficients  are  identified  by  solving  a 
linear  set  of  equations;  (b)  the  method  is  direct  (algorithm 
operates  on  the  data  directly,  without  the  need  to  estimate  the 
ACS)  ;  (c)  ARMA  (rather  than  just  MA  or  AR)  coefficients  are 
generated;  (d)  the  method  is  simple;  and  (e)  the  method  is 
implemented  efficiently  in  software. 

A  software  (MATLAB-based)  implementation  of  the  algorithm  was 
developed  and  exercised  using  simulated  multichannel  radar  data 
generated  with  the  SSC-developed  software  model  for  phased  array 
radar  in  an  airborne  surveillance  scenario  (Vol.  II  of  [Roman  and 
Davis,  1997]).  Simulation  results  obtained  in  the  study  validate 
the  modeling  capability  of  the  2D-LS-FD  algorithm,  and 
demonstrated  favorable  performance  as  a  whitening  filter  in 
comparison  with  the  conventional  optimal  joint-domain  method  using 
the  sample  covariance . 

In  Phase  II  a  MATLAB-based  software  package  will  be  developed 
to  test  STAP  algorithms  in  the  context  of  airborne  surveillance 
radar  arrays.  This  software  will  include  simulated  radar  data 
generation,  Monte-Carlo  analyses  capability,  and  extensive 
diagnostic  functions  (for  the  statistical  analysis  and  dissection 
of  intermediate  and  final  results) .  Classical  and  parametric  STAP 
algorithms  will  be  included  also.  This  software  will  be  applied 
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to  evaluate  the  MIBDA  methodology  based  on  2-D  hypothesis  filters, 
and  thus  establish  its  performance  in  relation  to  classical  STAP 
algorithms  and  to  STAP  methods  based  on  various  alternative 
parametric  algorithms  (Roman  and  Davis,  1993a,  1993b;  Michels, 
1991) .  The  software  will  be  generated  to  be  compatible  with  the 
structure  of  the  Rome  Laboratory  (RL)  STAP  (RLSTAP)  software 
analysis  and  simulation  package.  This  will  facilitate  integration 
of  the  software  into  RLSTAP  at  a  future  date. 

The  2D-LS-FD  method  will  be  the  baseline  2-D  algorithm  for 
Phase  II  due  to  its  attractive  .  features  and  the  performance 
established  in  Phase  I.  Alternative  formulations  of  the  method 
(such  as  image  noise  canceling)  will  be  evaluated  in  the  context 
of  the  RL  airborne  surveillance  phased  array  radar  application,  as 
well  as  the  to-be-selected  dual-use  application.  Model  stability 
is  an  important  issue  for  model-based  methods  in  general,  and  more 
so  in  the  case  of  2-D  models.  Other  candidate  algorithms  will  be 
implemented  in  MATLAB-based  software,  and  their  performance 
compared  with  that  of  the  baseline  algorithm. 
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APPENDIX  A. 


TWO-DIMENSIONAL  REPRESENTATIONS  FOR  ONE¬ 
DIMENSIONAL  MULTICHANNEL  STATE  SPACE  MODELS 


The  analytic  representation  of  the  sequence  transformation 
introduced  in  Section  2.1  (Equation  (2-1))  is  useful  in  the 
understanding  of  2-D  systems  and  of  the  differences/similarities 
between  1-D  and  2-D  systems  in  general.  Furthermore,  an  analytic 
representation  is  required  to  generate  the  2-D  scalar  transfer 
function  and  the  2-D  scalar  power  spectrum  of  a  2-D  system  related 
to  a  1-D  multichannel  system.  Therrien  (1981)  presented  such  a 
result  for  the  special  case  of  a  multichannel  AR  model,  and  it  is 
extended  herein  to  cover  state  space  models  in  innovations 
representation  form.  Availability  of  the  2-D  transfer  function 
and  power  spectrum  corresponding  to  a  1-D  multichannel  system 
allows  direct  comparison  of  the  frequency  domain  representation  of 
2-D  analytic  models  of  the  type  presented  in  Section  2-2  with 
state  variable  models  (SVMs)  of  the  type  employed  by  Roman  and 
Davis  (1993a;  1993b)  in  other  model-based  MIBDA  work. 

A .  1  Transformation  of  a  One-Dimensional  Vector  Sequence 

Into  a  Two-Dimensional  Scalar  Sequence 

Consider  the  1-D  multichannel  to  2-D  scalar  association  in 
Equation  (2-1),  repeated  herein  for  convenience: 

(A-l )  Xj.k(n)  =  x(n,k)  0<n<N-1;  0<k<J-1 

Recall  that  Equation  (A-l)  assumes  channel  J  is  the  temporal  and 
spatial  reference  for  the  array.  For  an  equally-spaced  linear 
array  with  array  spacing  d,  the  signal  at  channel  k  is  a  delayed 
(or  advanced)  version  of  the  signal  at  channel  k-1  .  With  channel  J 
as  the  temporal  and  spatial  reference,  it  follows  that 


where  fs  is  the  normalized  spatial  frequency  corresponding  to  an 
azimuth  angle  <|)  and  an  elevation  angle  0  (both  <|)  and  0  are  defined 
with  respect  to  the  array  boresight  for  a  side-looking  array) ,  and 
is  defined  as 

(A-3 )  fs=- — cos(0)sin(<(>) 

A,c 

where  X^  is  the  narrowband  radiation  wavelength.  Element  k  of  ej(fs) 
is  the  operator  representing  k-1  spatial  advances.  Vector  ej(fs)  is 
referred  to  as  the  spatial  frequency  vector  at  frequency  fs,  and 
can  be  expressed  in  terms  of  the  spatial  Z- transform  variable 
restricted  to  the  unit  circle  as 

1 

pj27rfs 

(A-4 )  0J(f.)=  ®  . 

ej2nfsK 


(A- 5)  zs  =  ej27tfs 
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Consequently,  the  spatial  frequency  vector  can  be  denoted  also  as 
ej(zs).  This  alternative  notation  is  more  appropriate  for  analyses 

involving  transfer  functions.  From  Equations  (A-l)  and  (A-2),  the 
association  between  a  1-D  multichannel  vector  sequence  and  a  2-D 
scalar  sequence  is  defined  as 

(A-6)  x(n)  =  x1(n)eJ(fs)  {x(n,k)  I  k  =  0,1 . J- 1} 

In  Relation  (A-6)  the  subscript  "1"  is  dropped  from  the  2-D 
sequence,  for  simplicity. 

A.  2  Innovations  Representation  State  Variable  Model 


The  innovations  representation  of  the  channel  output  process 
{x(n)}  is  a  linear,  shift -invariant,  stochastic  SVM  of  the  form 

(A-7a)  a(n+1)  =  Fa(n)  +  Ke(n)  n>0 

(A-7b)  x(n)  =  HHa(n)  +  e(n)  n  >  0 

(A-7c)  a(0)  =  0 

(A-7d)  E[a(n)aH(n)]  =  n(n)  =  n  n  *  0 


where  F  is  the  nsxns  system  matrix,  K  is  the  nsxJ  input 
distribution  matrix,  and  H  is  the  nsxJ  observation  matrix.  In  this 
model,  a(n)  is  the  ns-element  state  vector,  x(n)  is  the  J-element 
output  vector,  and  {e(n)}  is  a  zero-mean,  white,  Gaussian,  J-element 
vector  sequence  with  correlation  matrix  structure  given  as 


(A-8 )  E[e(n)8H(n-m)]  =  H  8(m) 


(A-9) 


m  =  0 
m  *  0 


V  n,  m 


48 


The  sequence  {s(n)}  is  referred  to  as  the  innovations,  and  matrix  K 
in  Equation  (A-7a)  is  the  Kalman  gain.  An  innovations 
representation  SVM  has  several  properties  and  characteristics  that 
simplify  model  identification  and  other  problems  (Anderson  and 
Moore,  1979;  Roman  and  Davis,  1993b) . 

Consider  now  an  instantaneous,  linear  transformation  on  the 
innovations  sequence  based  on  the  LDU  decomposition  of  the 
innovations  covariance  matrix,  £2.  That  is, 

(A-10 )  y(n)  =  THe(n) 

(A-ll )  th  =  d-1/2l1 

(A-12 )  n  =  LDLH  =  (T'1)Hr1 

In  the  LDU  decomposition  of  Equation  (A-12)  the  JxJ  matrix  D  is 
diagonal  and  real-valued,  and  the  JxJ  matrix  L  is  lower-triangular 
with  1's  along  the  main  diagonal.  The  covariance  matrix  of  the 
transformed  innovations  sequence  is  the  identity  matrix, 

(A-13 )  E[v(n)yH(n-m)]  =  I  8(m)  V  n,  m 

This  transformed  innovations  is  both  spatially  and  temporally 
uncorrelated  (Therrien,  1981) .  In  this  new  basis  for  the  input 
sequence,  SVM  (A-7)  becomes 

(A-14a)  a(n+1)  =  Fa(n)  +  KCv(n)  n>0 

(A-14b)  x(n)  =  HHa(n)  +  Cv(n)  n>0 

where  the  JxJ  matrix  C,  introduced  herein  for  notational 
simplicity,  is  defined  as 
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(A-15 ) 


C  =  (Th)'1  =(r1)H  =  LD1/2 


It  is  important  to  note  that  this  alternative  SVM  representation 
is  equivalent  to  SVM  (A-7)  from  an  input-output  perspective,  but 
SVM  (A-14 )  is  not  an  innovations  representation.  SVM  (A-14)  with 
{v(n)}  uncorrelated  and  normalized  as  in  Equation  (A-13)  is  referred 

to  herein  as  a  generalized  innovations  representation. 

The  JxJ  transfer  function  matrix  of  the  SVM  (A-14)  is  denoted 
herein  as  TQ|(zd),  where  Zd  is  the  temporal  (Doppler)  Z- transform 

variable  restricted  to  the  unit  circle, 

(A-16 )  zd  =  ej2ni(S 

and  fd  is  the  normalized  Doppler  frequency  shift.  In  terms  of  the 
model  parameters,  Tg,(  zd)  is  given  as 

( A- 17 )  TGI(2d)  =  TIR(2d)  C  =  [  HH  (Zjl  -  F)-1  K  + 1  ]  C  =  Hh  (zdl  -  F)’1  KC  +  C 

where  T)R(zd)  is  the  transfer  function  matrix  of  the  innovations 

representation,  SVM  (A-7) .  SVM  (A-14)  is  introduced  herein  to 
simplify  the  generation  of  2-D  frequency-domain  representations 
for  multichannel  SVMs,  as  shown  in  Section  A-4. 

The  square  matrix  factor  (zdl  -  F)'1  in  T|R(zd)  contributes  the 

multivariable  poles  of  the  SVM.  Specifically,  the  determinant  of 
matrix  (zdl  -  F)  is  a  polynomial  in  the  variable  Zd  whose  roots  are 

the  multivariable  system  poles.  In  contrast,  all  the  matrices  in 
the  transfer  function  expression  contribute  to  the  value  of  the 
multivariable  system  zeros  (the  multivariable  zero  definition 
preferred  herein  is  that  of  Davison  and  Wang  [1974,  1976]).  The 
dynamic  behavior  (in  the  temporal  domain)  of  the  SVM  is  determined 
by  the  multivariable  system  poles  and  zeros  jointly. 
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A .  3  Whitening  Filter  State  Variable  Model 


An  important  feature  of  the  innovations  representation  is 
that  SVM  (A-7)  is  causal  and  causally- invertible .  As  such,  the 
inverse  system  is  obtained  directly  from  SVM  (A-7) .  Specifically, 
the  inverse  system  for  SVM  (A-7)  is 


(A-18a) 

a(n+1 )  =  [F  -  KHh  ]  a(n)  +  Ke(n)  =  Aa(n)  +  Kx(n) 

n  >  0 

(A-18b) 

e(n)  =  -HHa(n)  +  x(n) 

n  >0 

(A-18c) 

a(0)  =  0 

(A-18d) 

E[a(n)aH(n)]  =  n(n)  =  n 

n  *  0 

where  all  vectors  and  matrices  are  defined  previously  except  the 
nsxns  matrix  A,  which  is  defined  implicitly  as 

(A-19)  A  =  F-KHh 

Notice  that  the  input  to  this  SVM  is  the  channel  output  sequence, 
and  the  output  is  the  temporally-uncorrelated  sequence  {e(n)} . 

Thus,  SVM  (A-18)  is  the  whitening  filter  corresponding  to  SVM  (A- 
7) .  An  important  relation  between  SVM  (A-7)  and  its  inverse  SVM 
(A-18)  is  that  the  multivariable  system  poles  and  zeros  of  SVM  (A- 
7)  are  the  multivariable  system  zeros  and  poles,  respectively,  of 
SVM  (A-18) .  Also,  SVM  (A-18)  is  an  innovations  representation 
(the  inverse  system  of  an  innovations  representation  is  itself  an 
innovations  representation) . 

Consider  again  the  spatially-whitening  instantaneous,  linear 
transformation  applied  to  the  innovations  sequence.  Equation  (A- 
10) .  Application  of  this  transformation  to  the  output  equation  of 
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the  whitening  filter  leads  to  the  following  generalized  whitening 
filter, 

( a-2 Oa )  a(n+1 )  =  [F  -  KHh  ]  a(n)  +  Ke(n)  =  Aa(n)  +  Kx(n)  n  >  0 

(A-20b)  v(n)  =  -  THHHa(n)  +  THx(n)  n>0 

As  before,  the  innovations  representation  property  is  lost  by  the 
introduction  of  the  linear  transformation  on  e(n). 

The  JxJ  transfer  function  matrix  of  the  SVM  (A-20)  is  denoted 
herein  as  TGW(zd),  where  Zd  is  the  temporal  (Doppler)  Z-transform 

variable  defined  previously,  Equation  (A-16).  In  terms  of  the 
model  parameters,  Tqw(  zd)  is  given  as 

(A-21)  TGW{zd)  =  Th  TWF(Zd)  =  TH[-HH(zdl  -  A)'1  K  + 1  ]  =  -ThHh  (zdl  -  A)'1  K  +  TH 

where  TWF(zd)  is  the  transfer  function  of  the  whitening  filter,  SVM 
(A-18 ) . 

A .  4  Two-Dimensional  Frequency-Domain  Representations  of  a 
Multichannel  State  Variable  Model 


Two-dimensional  frequency-domain  representations  are  derived 
next  for  each  of  two  SVMs :  the  generalized  innovations 
representation  (A-14),  and  the  generalized  whitening  filter  (A- 
20).  Consider  first  the  generalized  innovations  representation 
case.  The  input  to  the  generalized  innovations  representation  is 
the  generalized  innovations  sequence  {y(n)} .  A  relation  analogous 

to  (A-6)  can  be  established  also  for  the  generalized  innovations; 
namely, 

(A-22 )  v(n)  =v1(n)  fiJ(f8)  <^>  {v(n,k)  I  k  =  0,  1 . J-1} 
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Now  let  X(zd,zs)  and  N(zd,zs)  denote  the  2-D  Z-transform  of  {x(n,k)}  and 
{v(n,k)},  respectively.  Given  these  definitions  and  the  form  of  SVM 
(A-14),  the  scalar  2-D  transfer  function  TG|(zd,Zs)  of  the  1-D 

multichannel  system  (A-14)  is  obtained  as 
(A— 23 )  X(zd,zs)  =  Tel(zd,zs)  N(zd,zs) 

(A-24)  TGI(zd,zs)  =  -je»(zs)TGI(zd)eJ(zs)  =  le“(zs)T,R(zd)CeJ(zs) 

Figure  A-l  is  a  block  diagram  representation  for  these  relations. 
The  factor  J'1  in  Equation  (A-24)  is  required  to  preserve  the  power 
between  the  input  and  output.  Consider  the  case  where  TGl(zd)  is 
the  identity  matrix.  In  such  a  case  it  follows  that  {x(n,k)}  =  {v(n,k)}, 
provided  the  factor  J'1  is  included  (since  §j(zs)ej(zs)  =  J )  .  It  is 

possible  to  view  this  factor  in  a  different  way.  Specifically, 
the  spatial  frequency  vector  can  be  defined  to  be  of  unit  norm, 
which  requires  normalization  of  the  right  side  of  Equation  (A-4) 
by  a  factor  of  the  form  J'1/2. 


{v(n,k)} 


{x(n,k)} 


Figure  A-l.  Block  diagram  for  the  2-D  scalar  transfer  function  of 
a  2-D  system  related  to  a  1-D  multichannel  system  in  generalized 

innovations  representation  form. 
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The  temporal  frequency  content  of  the  2-D  transfer  function 
TG)(zd,Zs)  is  determined  by  the  temporal  dynamics  (multivariable 

poles  and  zeros)  of  the  generalized  innovations  transfer  function, 
tg,(  Zd).  In  turn,  the  spatial  frequency  content  of  the  2-D  transfer 
function  is  determined  by  the  spatial  frequency  vector  €>j(zs),  which 

inherently  models  a  delay,  but  does  not  include  spatial  dynamics . 
This  condition  may  limit  the  capability  of  1-D  vector  models  to 
represent  scenarios  where  the  coupling  between  channels  is  partial 
(where  the  noise- free  output  of  a  given  channel  is  not  a  delayed 
or  advanced  replica  of  the  noise-free  output  of  the  adjacent 
channels) .  Such  scenarios  include  cases  where  the  crab  angle  is 
non-zero  and/or  where  the  clutter  ridge  slope  is  different  from 
unity.  Further  analysis  of  the  modeling  capability  of  1-D  vector 
SVM's  is  required  to  assess  this  issue,  and  will  be  carried  out  in 
the  proposed  Phase  II . 

The  2-D  power  spectrum  of  the  linear,  shift-invariant  system 
in  Figure  A-l  is  the  square  of  the  magnitude  of  the  scalar 
transfer  function  evaluated  at  the  unit  circle  in  the  complex 
plane  (with  fd  and  fs  replacing  Zd  and  Zs  in  order  to  represent  this 

fact) ;  that  is, 

(A-25)  SGI(fd,fs)  =  TGI(fd,fs)T*,(fd,fs)  =  |TG|(fd,fs)|2 

-0.5  <  fd  <  0.5  and  -0.5  <  fs  <  0.5 

And  the  power  spectrum  of  the  2-D  channel  output  sequence,  denoted 
herein  as  Sxx(fd,fs),  is  determined  as  (Anderson  and  Moore,  1979) 

(A-26)  Sxx(fd,fs)  =  ^Gl^d’^s)  Svv(fd,fs)  =  |"^GI^d’^s)|  ^vv(^d’^s) 


-0.5<fd<0.5  and  -0.5<fs<0.5 
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where  Svv(fd,fs)  is  the  power  spectrum  of  the  2-D  generalized 
innovations  sequence,  {v(n,k)}.  Since  {v(n)}  is  temporally  and 
spatially  white  (Equation  (A-13)),  it  follows  that 

(A-27 )  E[v(n,k)v*(n-m,k-^)]  =  8(m,^)  V  n,  m,  k,  £ 


(A-28) 


8(m^) 


m  =  0  and  £  =  0 
m  *  0  or  £  *  0 


(A-29)  Svv(fd,fs)  =  1 


vfd)fs 


Then  Equation  (A-26)  becomes 

(A-30)  Sxx(fd,fs)  =  SGI(fd,fs)  =  |TGI(fd,fs)|2 

-0.5  <  fd  <  0.5  and  -0.5  <  fs  <  0.5 

which  is  the  desired  frequency-domain  representation  for  the 
generalized  innovations  SVM. 

Consider  now  the  generalized  whitening  filter  case.  The 
input  to  the  generalized  whitening  filter  is  the  channel  output 
sequence,  and  the  output  is  the  generalized  innovations  sequence. 

Based  on  Relations  (A-6)  and  (A-22)  and  on  the  form  of  SVM  (A-20), 
the  scalar  2-D  transfer  function  TGW(zd,Zs)  of  the  multichannel 

generalized  whitening  filter  is  obtained  as 
(A-31)  N(zd,zs)  =  TGW(zd,zs)  X(zd,zs) 

(A-32 )  TGW(zd,zs)  =  -1  ej(zs) TGW(zd)  ej(zs)  =  -j  e^(zs)  TH TWF(zd) eJ(zs) 
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with  TGW(zd)  as  defined  in  Equation  (A-21)  .  Transfer  function 
TGW(zd>zs)  i-s  presented  in  block  diagram  form  in  Figure  A-2  . 


Figure  A-2.  Block  diagram  for  the  2-D  scalar  transfer  function  of 
a  2-D  system  related  to  a  1-D  multichannel  system  in  generalized 

whitening  filter  form. 


The  2-D  power  spectrum  of  the  linear,  shift-invariant  system 
in  Figure  A-2  is  the  square  of  the  magnitude  of  the  scalar 
transfer  function  evaluated  at  the  unit  circle  in  the  complex 
plane  (with  fd  and  fs  replacing  zd  and  Zs  in  order  to  represent  this 

fact);  that  is, 

( A- 3  3  )  SQW(f d  >  fs )  =  ^GW^d  >  )  "*"Gw(f d  > *s )  =  l"*"Gw(f d  > )| 

-0.5  <  fd  <  0.5  and  -0.5<fs<0.5 

Equation  (A-33)  is  the  desired  frequency-domain  representation  for 
the  generalized  whitening  filter  SVM.  And  the  power  spectrum  of 
the  2-D  generalized  innovations  sequence  is  determined  as 
(Anderson  and  Moore,  1979) 
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(A-34 ) 


Sw(fd.fs)  =  SGW(fd.fa)Sxx(fd.fs)  =  lTGw(fd.U|2  Sxx(fd,fs) 

-0.5  <  fd  <  0.5  and  -0.5  <  fs  <  0.5 

with  Svv(fd,fs),  the  power  spectrum  of  the  2-D  generalized  innovations 
sequence,  of  the  form  in  Equation  (A-29) . 
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OF 
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Mission.  The  mission  of  Rome  Laboratory  is  to  advance  the  science  and 
technologies  of  command,  control,  communications  and  intelligence  and  to 
transition  them  into  systems  to  meet  customer  needs.  To  achieve  this, 
Rome  Lab: 


a.  Conducts  vigorous  research,  development  and  test  programs  in  all 
applicable  technologies; 

b.  Transitions  technology  to  current  and  future  systems  to  improve 
operational  capability,  readiness,  and  supportability; 

c.  Provides  a  full  range  of  technical  support  to  Air  Force  Material 
Command  product  centers  and  other  Air  Force  organizations; 

d.  Promotes  transfer  of  technology  to  the  private  sector; 

e.  Maintains  leading  edge  technological  expertise  in  the  areas  of 
surveillance,  communications,  command  and  control,  intelligence, 
reliability  science,  electro-magnetic  technology,  photonics,  signal 
processing,  and  computational  science. 

The  thrust  areas  of  technical  competence  include:  Surveillance, 
Communications,  Command  and  Control,  Intelligence,  Signal  Processing, 
Computer  Science  and  Technology,  Electromagnetic  Technology, 
Photonics  and  Reliability  Sciences. 
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